Gaussian approximation to single particle correlations at and 
below the picosecond scale for Lennard-Jones and 

nanoparticle fluids 

R. van Zon 1 ' 3 , S. S. Ashwin 2 ' 3 and E. G. D. Cohen 3 

00 ■ 
O , 

O ! 1 Chemical Physics Theory Group, Department of Chemistry, 

(N 



oo 



University of Toronto, 80 St. George Street, Toronto, Ontario M5S 3H6, Canada 

2 Department of Chemistry, University of Saskatchewan, 
110 Science Place, Saskatoon, Saskatchewan S7N 5C9, Canada and 

3 The Rockefeller University, 1230 York Avenue, 
New York, New York, 10065-6399, USA 
(Dated: 28 March 2008) 

To describe short-time (picosecond) and small-scale (nanometre) transport in flu- 
ids, a Green's function approach was recently developed. This approach relies on 
an expansion of the distribution of single particle displacements around a Gaussian 
function, yielding an infinite series of correction terms. Applying a recent theorem 
[Van Zon and Cohen, J. Stat. Phys. 123, 1-37 (2006)] shows that for sufficiently 
small times the terms in this series become successively smaller, so that truncat- 
^ \ ing the series near or at the Gaussian level might provide a good approximation. 

£h ' In the present paper, we derive a theoretical estimate for the time scale at which 
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truncating the series at or near the Gaussian level could be supposed to be accu- 
rate for equilibrium nanoscale systems. In order to numerically estimate this time 
scale, the coefficients for the first few terms in the series are determined in computer 
simulations for a Lennard-Jones fluid, an isotopic Lennard-Jones mixture and a sus- 
qq | pension of a Lennard-Jones-based model of nanoparticles in a Lennard-Jones fluid. 

The results suggest that for Lennard-Jones fluids an expansion around a Gaussian 
is accurate at time scales up to a picosecond, while for nanoparticles in suspension 
(a nanofluid), the characteristic time scale up to which the Gaussian is accurate 
OO ' becomes of the order of five to ten picoseconds. 

o 

PACS numbers: 05.20.-y, 02.30.Mv, 02.60.Cb, 61.20.Ja, 05.60.Cd 
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I. INTRODUCTION 

Small clusters of particles suspended in a fluid occur in many forms, from nanoparticles 
[J 0, 01, quantum dotsj^] and colloidal suspensions jl] to biomolecules such as globular pro- 
teins ja, [l]. Such nanoclusters have a variety of applications, from material coatings to d rug 



delivery by hollow clusters. Both the individual behaviour of nanosized particles [9|, [1CJ, III 
as well as their collective behaviour, such as the increased heat conductance in dilute sus- 
pensions of nanoparticles (so-called nanofluids) []]] , have received considerable attention^. 

For the purpose of studying small length scale and short time classical transport phenom 
ena which occur in nanosystems, a Green's function approach was introduced by Kincaid[F" 



This approach has the promise of being able, in principle, to describe transport phenomena 



2 



on all time and length scales, unlike hydrodynamics. The main idea of the theory is to 
describe the evolution of fluid properties such as its energy, momentum and number density 
in terms of Green's functions. The application of these Green's functions to nanosystems 
and systems where time scales at picoseconds or less are important, has been an area of 



some interest [131. Il4j . 1 151 . Il6| . In these cases, the Green's functions were expanded around 
a Gaussian distribution plus an infinite series of corrections, a finite truncation of which 
yielded excellent agreement with simulations. Even just the Gaussian itself was found to 
be a reasonable approximation to the Green's functions. An explanation for this could be 
that the series of corrections has fast convergence, but at that point, it was not known why 
that this could be the case. Since the Gaussian description is much simpler than the full 
Green's function, one would like to know when fast convergence occurs and when taking the 
Gaussian approximation suffices. A preliminary answer to this question was found in Ref.ll7l. 
namely, that for the motion of a single particle in an equilibrium pure Lennard- Jones (LJ) 
fluid, the Gaussian approximation can be used up to time scales of the order of a picosecond. 

One of the applications of the Green's function approach is mass transport in liquids 
and liquid mixtures. For that case, the Green's functions are essentially the probability 
distribution functions of displacements (in a time t) of single particles of the different com- 
ponents [lq |. Thus it is not too surprising that the Green's functions can be expressed in 
terms of the cumulants of this distribution. These cumulants measure the correlations of the 
displacement of a single particle, in particular, they measure the departure of the correla- 
tions from Gaussian behaviour. As will be discussed in more detail below, a recent theorem 
regarding these cumulants implies that when the Green's functions are expanded around 
a Gaussian distribution, the correction terms to the Gaussian term are proportional to in- 
creasing powers of t for short (initial) times t 18[ . Analytic expressions for the coefficients 
in front of the powers of t were also derived in Ref.0. The values of the first two numerical 
coefficients are here of particular interest, because they can be computed numerically and, 
as show in Sec. lVBI can then be used to find estimates of the physical time scales below 
which the expansion of the Green's function around the Gaussian term yields useful results, 
as appeared to be the case in Refs.El, EH, EH EEl Numerical values for these coefficients will 
be presented in this paper for various equilibrium LJ-based systems, including nanoparticles 
in a suspension of LJ particles. We will present the resulting orders of magnitude of the 
relevant time scales on which the first few terms in the series decrease. Non-equilibrium 
systems will be studied in future work. 



II. SYSTEMS 

Three systems were studied, namely a pure LJ fluid, an isotopic binary mixture of LJ 
particles (in which context the study of short time displacements arose [3]), and a suspension 
of nanoparticles in a LJ fluid. 

In the isotopic binary LJ mixtures, there are Na particles of mass tua and Nb particles 
of mass rriB in a box of size L 3 , such that the number density is p = (Na + Nb) / L 3 . For the 
pure LJ fluid, one sets Nb = 0. The positions and velocities of the particles will be denoted 
by r\i and v\i, respectively, where A = A or B and i is a particle index, which runs from 1 
to Na if A = A and from 1 to Nb if A = B. By definition, in an isotopic mixture all pair 
interaction potentials are the same for all components, but their masses are different. The 
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inter-atomic potential between the particles is the LJ potential 

V AA (r) = V AB (r) = V BB (r)=4e 



<7\ 12 

r ■ 



(1) 



where r is the distance between two particles, and a and e are the same for all pairs of 
particles. 

All quantities reported are in LJ units: length in units of a, temperature in units of e/k B , 
density (p) in units of a" 3 and time in units of tlj = [a 2 m A j t) x l 2 ^ where m A is the mass of 
an A-particle. In other words, we will use units in which a = 1, k B = 1, e = 1, and m A = 1. 
Although these are arbitrary units, to understand the physical consequences of our results, 
we use the LJ parameters of Argon as a reference. In that case, one unit of time corresponds 



to t B j = 2.16 x 10 12 seconds, while one unit of length corresponds to o = 0.34 nm 19|, |20 



As mentioned above, apart from the pure L J fluid and the isotopic binary L J fluid mixture, 
a third system which will be studied, namely, a suspension of nanosized particles in a fluid, 
often called a nanofluid. One can obtain this system from the binary isotopic LJ fluid 
mixture by changing the B particles to much larger, nanosized particles while the A particles 
remain regular LJ particles, and changing the potentials V AB and V BB in the following way. 
Each nanoparticle is represented as a spherical cluster of radius R with a smoothed uniform 
distribution of M LJ particles as proposed in Refs.0 and 21. Since we are only after typical 
time scales for which the expansion presented in Sec. lIIII below is valid, we restrict ourselves 
here to this simple nanoparticle model. For simplicity, we therefore take the strength of the 
LJ potential between the constituent LJ particles of the nanoparticles and the fluid particles 
to be the same, and the mass of the constituent LJ particles of the nanoparticle is also taken 
to be equal to that of the fluid particles. R will range from 1 to 6 in LJ units, i.e. from 0.34 
nm to 2 nm (which is a typical size of a quantum dotj^j), while M will be chosen such that 
for R = 0, the nanoparticle reduces to a single LJ particle (M = 1) while for large R the 
density of L J particles within the nanoparticle approaches one. This can be accomplished by 
choosing M to be 1 + -R 3 , leading to a maximum mass ratio of 217 between the nanoparticles 
and the fluid LJ particles. One can show that the result of integrating the LJ potentials 
corresponding to all the points in the spherical nanoparticle is that a nanoparticle interacts 
with a fluid LJ particle through the potential^, 2l| 



V AB {r) = AM 
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where r is the distance between the centre of the nano particle and the LJ fluid particle, 
while the interaction potential between two nanoparticles is given by(2l| 



V BB (r) = 4M 2 
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where r is the distance between the centres of the nanoparticles. Note that because of the 
much larger size of the nanoparticles, far fewer will fit into a system of given volume than 
B particles fit in an isotopic LJ mixture of only LJ particles. 
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The systems studied in this paper are all in canonical equilibrium, i.e., their distribution 
function p eq (X) in phase space (T = {r A j,v Ai }) is given by: 

p eq {T) = e'^/Z, (4) 

where Z = j exp[— H(r)/T]dT is the partition function, T is the temperature, and H is the 
Hamiltonian which is of the form 



TV 



H{T)= £ y;T^ + u, ( 5) 

X=A,B j=l 

where U is a sum of pair potentials: 

U = E E E E l^(|r A ,-r w -|), (6) 

X=A,B i=l ti=A,B j=l 

where the prime excludes equal particles (i.e., A = fi and i = j) and the V X/1 are of the form 
given in Eqs. (Op)-® above. Finally, we remark that the equations of motions are given by 

1 dU 

r\i = v A ,; v Ai = — . (7) 

m x dr Ai 

III. GREEN'S FUNCTIONS AND CUMULANTS 

We will now briefly review the Green's functions approach and its connection with the 
distribution of single particle displacements. For mass transport processes, the number 
density n A (r,t) of a specific component A at position r at time t can be written as 

n x (r, t) = J dr' G x (r - r', r', t)n x (r', 0), (8) 

where G\(r,r',t) is the Green's function for component A (A or B for a binary mixture), 
which is defined aspj], lfl 



(d[r' - r Ai (0)J) is 

where r X i(t) is the position of the ith particle of component A at time t and the average 
()i s is over a (possibly non-equilibrium) initial state ("is"), which has to be specified for 
the particular problem that one wants to study. The Green's function G A (r,r',t) can be 
interpreted as the probability that particle i of component A was displaced over r in a time 
t given that it started at r'. Note that the Green's functions do not depend on i because 
particles of the same kind are indistinguishable. 

Although the Green's function approach is aimed primarily at non-equilibrium systems, 
we will restrict ourselves here only to equilibrium systems, because the time scales for the 
validity of the expansion to be presented below are expected to be similar in equilibrium 
and not-too-far-from-equilibrium systems, and the equilibrium system is much easier to deal 
with from a numerical point of view. In the equilibrium case, the Green's functions become 
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independent of r' because the system is homogeneous and are then identical to the Van Hove 
self-correlation functions G A (r, t) (with A a component) defined as [221 



1 Nx 

G X s (r, t) = W J2 + r ><(°) " '*(*)]> > ( 10 ) 

i=l 

where the subscript s refers to G A being a self-correlation function of a single particle. The 
average ( ) is here taken over the canonical equilibrium ensemble p eq given in Eq. (Tj0). To 
see that Eq. ( fTOl) is the equilibrium variant of Eq. ([9]), note that each term on the right-hand 
side of Eq. (1101) gives the same contribution to the sum due to the indistinguishability of 
particles of the same component. Thus one can also write 

G^(r,t) = (5[r + r A1 (0)-r A1 (t)]), (11) 

where particle 1 of component A is used as a representative particle of that component. The 
expression for the Van Hove self-correlation function in Eq. ( 11 II) coincides with that for the 
Green's function in Eq. ([9]) in cases where the Green's functions have no r' dependence, 
i.e., in equilibrium. Note that like the Green's function, the Van Hove self-correlation 
function G A (r, i) can therefore be interpreted as the probability that a single fluid particle 
of component A has experienced a displacement r in a time t. 

The Fourier transform of the Van Hove self-correlation function is the self-scattering 
function F A (k, t) j22|, which is given by: 

F*(k,t) = ( e ifck -I r ^- r ^°)]) = ( e ikAx ^) . (12) 

Here k = kh is a wavevector with length k along the unit vector k and 

Ax A1 (t) = k-[r A1 (t)-r A1 (0)] (13) 
denotes the displacement of particle 1 of component A along the direction k at a time t. The 



self-scattering functions can be measured by incoherent neutron scattering experiments 23 

According to elementary probability theory j24j one can interpret log F^(k,t) as the cu- 
mulant generating function of Ax\i{t), where Ax\i(t) is considered to be a random variable, 
so that F*(k,t) can be written in the following form: 



F s \k,t) =exp 



oo A 
n=l 



(14) 



Here k a is called the nth cumulant of the displacement Ax\i(t). The behaviour of these 
cumulants as a function of time has been investigated in the context of incoherent neutron 



scattering by Schofield 25] and Sears 26]. They showed that for equilibrium systems, the 
cumulants (k„ for n = 2, 4, 6) have the following behaviour at small times: ki ~ 0(t 2 ), 
«4 ~ 0(t 8 ) and k,q ~ 0(t 12 ), while the odd cumulants vanish in equilibrium. This behaviour 
suggested a generalization, which has recently been obtained for a certain class of physical 
systems as a Theorem [181 ] . For a class of classical systems which includes systems with 
smooth potentials 1 in canonical equilibrium, it was shown that the K^(t) have the following 



1 The LJ potential is not truly smooth because it diverges at r = 0. However, in equilibrium, this point 
has a vanishingly small probability, so that the LJ potential may be treated as effectively smooth. 
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form: 



c x t n + 0(t n+1 ) 
c x t 2n + 0{t 2n+l ) 



for n < 3 
for n > 3. 



(15) 



where c„ are coefficients independent of t. We see from Eq. (1141) that for sufficiently small 
wavevectors k, F x (k,t) ~ exp[— k x k 2 /2]. Since F x (k,t) is then approximately Gaussian in 
k, we would expect that its inverse Fourier transform, the Van Hove self-correlation function 
G x (r,t), is also approximately Gaussian in r. The corrections to the Gaussian behaviour 
of F x (k,t) are given by the terms in the series in Eq. (1141) with n > 2. Taking the inverse 
Fourier transform of Eq. (|14p . one can show that the Van Hove self-correlation function is of 
the form of a Gaussian plus corrections 18]: 



G x (r,t) 



exp( 



-w 



1 + 



k x H 4 {w) k x H 6 (w) 



A\A[k x } 2 
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6!8[41 3 
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(16) 



Here H n is the nth Hermite polynomial, and w = r j \J2k x a dimensionless length. Substi- 
tuting Eq. (1151) in Eq. f)16H . the Van Hove self-correlation function can be expressed as a time 
series of the form: 



G x (r,t) 



exp( 



-w 



c x m\t & 
5760T 3 



H 6 (w) + 



(17) 



where we used that in equilibrium c x = (v^) = T/m x . 

There are a few systems for which all the c x for n > 2 are zero, leading to Gaussian 
Van Hove self-correlation functions. These systems are the ideal gas and systems with only 
harmonic forces, whose equations of motion are linear. For nonlinear systems, however, the 
right-hand side of Eq. (1171) is a series in increasing even powers of t. It is natural to expect 
that for a small enough t, the successive terms in these series should rapidly decrease. This 
would mean that the series converges and that one could use a finite number of terms, or 
even just the Gaussian, as a good approximation to the whole series. Applying the general 
rule that a series ^^L a n converges if lim^oo \a n+ i/a n \ < 1 to the series in Eq. (fTTI) . where 
a n oc C2 n t 2n , it follows that the time scale below which the decrease in the terms occurs 
depends critically on the coefficients c x n , or in particular on ratios of successive c x n as n 
approaches infinity. Infinitely large values of n are, of course, beyond the reach of numerical 
computation but to get an estimate for the time scales, we numerically evaluated c^'s for 
the LJ liquid for finite n up to n = 3 and the corresponding time scales for the decrease in 
the terms of the series. 



IV. TIME SCALES 

As explained above, to numerically estimate the time scales up to which the series ex- 
pansion of the Van Hove self-correlation functions G x (with A = A or B) in Eq. (j!7p may 
converge or at least be useful, we are interested in the first few terms of the series. The 
terms in Eq. (JT7J) which are of importance are then the coefficients c\ and c x . Expressions 
for these coefficients are derived in Sec. El while in Sec. lVIl the results of their numerical 
evaluation in simulations are presented. 

For sufficiently small times t, every successive term in the series in Eq. (II 7p would ap- 
proach zero more rapidly than the previous term because of a larger power of t associated 
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with it. This gives us a simple relation to check when we could expect the terms in the 
series to decrease. The first estimate of a time scale, to be denoted by t g , follows from the 
criterion that for t = t g , the first term in the brackets in Eq. (1171) . i.e. 1, is of the same order 
of magnitude as the next term, i.e. c^m^H^iw) j '(96T 2 ). To find the order of magnitude of 
the latter, we need an order of magnitude estimate for H^w), which we find as follows. The 
prefactor e~ w in Eq. ffTTl) suggests that w = 0(1), since otherwise G A would be extremely 
small. The Hermite polynomial H^(w) contains no physical parameters, only numerical fac- 
tors which are also of 0(1), so we conclude that H^w) = 0(1). The second term in Eq. f|T7|) 
is therefore of the order of the first term at t = t g with c^ml^] 4 / (96T 2 ) = 0(1), yielding 



This Tq expresses on what time scale a Gaussian approximation to G A will break down, 
while for time scales somewhat less than to t g , the Gaussian distribution could be supposed 
to be a good approximation. 

The next simplest estimate of a time scale, to be denoted by r A , is determined by the 
time t — t a when the second and third terms in the square brackets in Eq. flTTj) become 
comparable, i.e., when: 



c\m\t 
96T 2 
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HAw) 



5760T 3 " 



(19) 



which, using the same argument as above Eq. f fTHl) to show that typical values of H^{w) and 
Hq(w) are 0(1), leads to 

1/2 



60 
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1 1 




c 6 





T 

— ■ (20) 



This t£ also defines a time scale below which the subsequent terms in the series in Eq. (j!7p 
should decrease in magnitude. Thus, for time scales sufficiently less than t£, the Cg term 
can be neglected compared to the c\ term in Eq. ffTTl) . but for time scales larger than r^, the 
Cg term certainly needs to be taken into account. 

One could in principle get additional time scale estimates by including higher order 
terms in Eq. (|17|) and comparing the nth with the n + 1st term. Note that then t g is equal 
to t± and is equal to r^ 1 , respectively. If the limit r A = lim^oo-r^ exists, the series in 
Eq. (I17p converges for all t < r A . In simulations, we cannot take this limit, but we will 
see that r G and r A have similar orders of magnitude, suggesting that t g and r A might be 
reasonable estimates of the actual time scale of convergence of Eq. ffTTl) . 

V. EXPRESSIONS FOR THE COEFFICIENTS c A AND c A 
A. General expressions 

We first discuss the analytical expressions for the coefficients c A in terms of the so- 



called multivariate cumulants based on Ref.ll8l. The general relation between moments and 



cumulants is given infA] For short times, the K^(t) have the form given by Eq. ffT5|) . where 
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for n > 3 the scaling coefficients are given by[fj 

n n . 

' -E-E r T - +1 ,:^»„, (( y i" l| ;...;yfc7 1 )). (2D 



ni=0 n n+ i=0 -l J-7 



n::iKK7!) ? 



E;iiX=n 
E"S7"7=2n 

Here, ((E^ 1 '; • • • ; Exn+i 1 ')) * s a notation introduced in Ref.0 for a multivariate cumulant, 
which is a multivariate moment with all possible factorizations subtracted. In this notation, 
quantities separated by semicolons are treated as separate random variables and if a quan- 
tity has a superscript within square brackets, it denotes the number of repetitions of that 
particular quantity, e.g., ((Y^f)} = ((E A1 ; E A1 ; E A1 )) (see[S}. Furthermore, E A7 is defined as 



Y) 



A 7 



<Pkx xl (t) 



dP 



(22) 

i=0 



with Ax\\(t) defined in Eq. (fTBl . Note that we deviate here from the notation in Ref.lla. 
where the cumulants were expressed in terms of X x -y = Evy/7! instead of in terms of E A7 . 

By writing out the sums in Eq. (121 p for n = 4 and n = 6, one finds the following 
expressions for c\ and Cg: 

= M Y S^)) + iriM 1 )) + H Y S; n 4 » + ±«e ai; e| ] ; e A3 » + ^rS)) (23) 
4 = 8io(( y i?; y Ar)) + M*i 41 ; Y * ^» + M Y lh Y ^ y*)) + M Y ^ Y 3)) 

+ l«yg;Y$\ YS)) + i6 {{Y\i ; Eg; E A3 )) + ^((E^; E A2 ; E A3 ; E A4 )). (24) 

To evaluate these expressions, we need the explicit expressions for the E A7 . Since the E A7 
are simply the 7th derivative of Ax A1 , they can be found by straightforward differentiation 
(cf. Eqs. ([3) and (TI3"]) ). The resulting expressions are polynomials in the velocities of the 
p articles [181]. Below, it will turn out that only the highest power of the velocities in the 
expression of each E A7 leads to a non-zero contribution to c\ and Cg. It suffices therefore to 
write only the highest powers in the velocities for the E A7 , i.e., 

Y\\ = v xlx (25) 



E 



A3 



m\ dx xl 

E ft a — ■ v ^ 27 



E 



A4 



E 



A5 



jEE fe ft r :v^ + 0(^) (28) 

E E E T r : v^v* + <V) (29) 
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£tj k£ pn rp n rip 

where each sum over two indices denotes a sum over the components A and B for the Greek 
index and a sum over the particles of that component for the Latin index, while 0(v n ) 
represents terms which are a polynomial of order n in the velocities. 



B. Simplifications for equilibrium systems 

In equilibrium, the velocities are independent Gaussian distributed variables with zero 
mean (cf. Eqs. fl2J) and (jSD), which allows some simplifications in the expressions for c\ and 
Cg in Eqs. (j23p and (1241) . respectively. These simplification will not only lead to shorter 
expressions but will also reduce the number of quantities inside each cumulant, i.e., it will 
reduce the order of the cumulants. This is numerically advantageous since higher order 
cumulants tend to require more statistics to keep the error small. 

The first simplification is that, given the Gaussian nature of the velocities, Theorem A 



of Ref.ll8l can be applied to show that the terms denoted by 0(v n ) in Eqs. ff28]) - fl3T|) do not 
contribute to the right-hand side of Eqs. (T2"3"j) and flUJ), because they contribute cumulants 
which contain fewer powers of the velocity than the number of velocity factors Y\\ = v x \i in 
the cumulants, and according to Theorem A, such cumulants are zero (see the Appendix in 



Ref.ll8l for details). On the other hand, the first terms on the right-hand sides of Eqs. (128]) - 
( |3T|) contain just enough powers of the velocities to match the number of factors of Y\i = v x \i 
in the cumulants in Eqs. (123]) and (1241 so that Theorem A does not apply and they might 
yield a non-zero result. Thus only these terms in Eqs. (l28]) - (l3~Tj) need to be taken into 
account. 

The next simplification involves the average over the velocities, which can be taken sep- 
arately from the average over the positions because of the factored form of the canonical 
equilibrium distribution given in Eq. (j4"]). Thus, canonical averages can be taken in two steps: 
first an average over velocities and then an average over positions. To apply this two-step 
process to cumulants, one needs to relate the cumulants to averages. Using Eq. (1A3|) . the 
cumulants on the right-hand sides of Eqs. (123]) and (124"]) can be written in terms of moments 
which are simply averages of products of factors of Y\ T For velocity averages of products 
of independent Gaussian distributed velocities with zero mean, we can use Wick's theorem 
which states that the average can be obtained by pairing the velocities in all possible ways 
and then taking the average for each pair separately. Note that the average of two velocities 



v mii an d v^ 2 i 2 is 



T 

vViii^/xaia/u — ^^i/j,2 ^u*2*> (32) 



where the subscript v of the brackets indicates that only the average over velocities is per- 
formed. Afterwards, the average over positions, denoted by () r , still needs to be performed 
to obtain the full average. 

The straightforward method of writing the cumulants out in terms of moments intro- 
duces a lot of subtractions terms, which can be largely avoided by formulating a similar 
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Wick's rule for cumulants. However, the two-step nature of the averaging process, in- 
volving velocity as well as position averages, is a complicating factor here. Forgetting for 
the moment about the position average, for Gaussian distributed velocities, cumulants can 
be computed similarly as averages, i.e. using Eq. (1521) . with the distinction that there be 
only "connected contributions", in the sense that the pairing of velocities be such that 
all expressions in the cumulant are connected to each other. To give an example, for the 
cumulant ((viVj-,VkVi)) v , the term (viVj) v (vkVi) v does not connect the expressions V{Vj and 
VkVi, and therefore does not contribute, while the terms (viVk) v (vjVi) v and (viVi) v (vjVk) v 
do connect the two, so that ((viVj]VkVi)) v = (viVk) v (vjVi) v + (viVi) v (vjVk) v . However, when 
averaging with p eq in Eq. (J4J), there is a second, non-Gaussian, average, namely, over the 
positions. As a consequence, although a term like ( •£ ¥ ViVj) v ( £ V VkVi) v may seem dis- 



connected and therefore not to contribute to the cumulant (( q t X v i v j '■> aT k dr v k v it)i ^ ne sec_ 
ond average over positions will, as it were, reconnect the parts. One can show such 
seemingly disconnected expressions (as far as the velocities are concerned) still yield a 
contribution to the cumulant which is equal to the position-cumulant of the factors, i.e. 
(((aS^^' (aS^^)"))*' = S^))r( v i v j)v(vkVi) v , where a subscript r denotes a 

cumulant over the positions only. 

With these rules on how to compute cumulants, we now return to the expressions for 
C4 and Cg in Eqs. ( 1231) and ( j24j) , respectively. One easily checks that to get connected 
contributions, all the factors Yxi = v x \i in the cumulants in Eqs. ( 1231) and ( |24"1) must be 
paired with velocities in the other Y\ T If ri\ is the number of factors of Y\i m a cumulant, 
this introduces a factor n\\ due to the number of ways one can pair two sets of n\ velocities. 
Furthermore, because of the Kronecker delta's in Eq. (1321) . all summations from Eqs. (J2j 
( l3Tj) can easily be performed, and one finds 
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Here the same notation has been used as explained below Eq. ff2T]) and in [A] 
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The above expressions can still be further simplified for systems in canonical equilibrium, 



using the following identity due to Yvon[28l. 129 
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I dB 
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(35) 



for any function B of the position of the particles, as can be proved by partial integration. 
While we will not present the lengthy details here, this identity can be used to find linear 
relations between the expressions on the right-hand sides of Eqs. ( 13 6 j) and (1371) . which allow 
us to rewrite the expressions for c\ and Cg in a variety of ways. Among those, we choose 
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(36) 
(37) 



These equations require at most second and third order cumulants, respectively, which is 
advantageous since numerically higher order cumulants tend to produce larger statistical 
errors. They agree with the expressions found by Sears for a one-component fluid (26[. Note 
that in the special case of a harmonic potential, derivatives higher than the second vanish, 
so that then for and Cg only the last terms in Eqs. (1361) and (1371) . respectively, remain, 
which only involve the cumulants of the second derivative of the potential. Since the second 
derivative is constant for a harmonic potential, these cumulants are zero as well, so that the 
coefficients c\ and Cg are zero, as expected for a linear system. 

With this background, next, we will present the results of the numerical evaluation of the 
coefficients c\ and Cg for a number of equilibrium systems by means of molecular dynamics 
simulations, in order to estimate the time scales Tq and which indicate where one could 
suppose that the first term alone (i.e. the leading Gaussian) or the first few terms (i.e. the 
Gaussian plus corrections) of the series in Eq. (jPTl) can be used as a good approximation to 
the full Van Hove self-correlation function. 



VI. SIMULATION RESULTS 
A. Single component Lennard-Jones fluid 

In this section, we present the numerical result for C4 and cq (cf. Eqs. ( 1361) and (1371) ) and 
the resulting time scales tq and r* (cf. Eqs. (fTSl and (|20|) ) for a single component fluid of 
N = Na LJ particles with periodic boundary conditions in a box of linear size L = 5 (in LJ 
units). Note that we have omitted the component-superscript A here because there is only 
one component. The results were obtained from molecular dynamics (MD) simulations, for 
which the initial conditions were drawn from the canonical distribution by employing an 
isokinetic Gaussian thermostat [271] during the equilibration stage, while the runs themselves 
were done at constant volume and energy. In the simulation, a potential cutoff of r r = 2.5a 
was used and the equations of motion were integrated using the Verlet algorithm 1^] with 
a time step of 2 femtoseconds. 

Since tq and r* will depend on temperature and density, it is of interest to study the 
dependence of C4 and cq as a function of these two parameters. We studied the temperature 
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FIG. 1: The coefficients C4 (on the left) and c% (on the right) as a function of temperature T 
for an equilibrium single component LJ fluid with density p = 0.8. These results are from a MD 
simulation with N = 100 particles, with periodic boundary conditions. All quantities are in the 
LJ units defined in Sec.HIl 



dependence by keeping N and p fixed to 100 and 0.8, respectively, while temperature values 
ranging from 1 to 3 were used. For each of these parameter values, data were accumulated 
once equilibrium had been attained in the simulation and collected every 2 ps in a 8 ps 
long run, yielding five points per run. This was repeated for 2000 different initial conditions 
(yielding 10,000 points per temperature) for each temperature value and the results for C4 
and cq were averaged over these 2000 runs. To decrease the statistical errors even further, 
we averaged over all particles of the same kind (i.e. replacing the index 1 in Eqs. (1361) and 
( 13"T|) by any index i and averaging the results) as well as over the three directions of space 
(i.e. replacing x by y and z in Eqs. ( I3"rj|) and (157)1 and averaging). 

The resulting behaviour of C4 and c§ as a function of temperature is shown in Fig.[TJ The 
data for C4 in the left panel of Fig.[T] are consistent with the preliminary data that were 
presented in Ref.ll7l. Note that in Fig.[T], the absolute value of the coefficient cq has been 
plotted. The reason is that the values of cq that are found in the simulations are always 
negative. In Fig. (2J we plotted the resulting time scales tq and r* (cf. Eqs. ( !T8~]) and (1201 ) 
as a function of temperature. We see that by increasing the temperature, we moderately 
decrease these time scales from roughly 2 ps to 1 ps, which are the estimates for the time 
scales up to which the series in Eq. (fT7j) could be supposed to give an accurate approximation 

The density dependence of c 4 and c 6 was also investigated using the same setup, but 
keeping the temperature fixed at T = 1.0, while the density ranged from p = 0.5 to p = 1.0. 
The resulting time scales tq and function of density are plotted in Fig.[3j While 

both timescales remain on the order of one or two picoseconds under changes of the density, 
we see that the two time scales tq and r* behave quite differently; whereas the time scale 
tq decreases moderately with increasing density, indicating that the first correction term 
in Eq. (fTTj) becomes important somewhat sooner for higher than for lower densities, the 
time scale r* is virtually constant as a function of density and bigger than tq, indicating 
that the second correction term in Eq. (1171) becomes important at a slightly larger time scale. 
However, the order of magnitude of these two time scales is so similar (i.e. both of picosecond 
order) that such a distinction does not appear to be significant. 
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FIG. 2: The critical time scales tq and at which the series in Eq. (|17p for the Van Hove self- 
correlation function of an equilibrium single component fluid could be supposed to be practicable 
(cf. Sec.lFVj below Eqs. (fl8|) and (|20|) ) as a function of temperature T for a density p = 0.8. Note 
that the physical time scales in picoseconds can be calculated by multiplying both r* and tq by 
the LJ unit time tlj = 2.16 ps. 



B. Isotopic Lennard- Jones Binary Mixture 



Our investigation into the cumulants originated in the study of mass transport in binary 



isotopic mixtures at short time scales [16 



and hence we are interested in the time scales 
in binary isotopic mixtures as well. From the expressions for the time scales in 
Eqs. ( ITS]) and (1201) as well as for the coefficients c% and Cg in Eqs. ( 1331) and (jMj) . respectively, 



Tg and r. 



one can readily deduce that c\ oc m^" 4 and Cg oc m^ b . Using this in Eqs. (|T8l) and (1201) . 
one sees that the time scales r G and simply scale as the square root of the mass. The 
remaining parts of the coefficients only involve the potential, which in an isotopic mixture 
is the same as for a pure LJ system. Therefore, no new simulations are needed for this case; 
the time scales are those of the pure LJ system, multiplied by the square root of the mass 
ratio of the components and the original LJ particles, i.e.: 



1 G 



mx 
m 

mx 
m 



(38) 
(39) 



where m is the mass of the particles in a single component LJ fluid. 

Since in Nature, there are no isotopes with large mass ratios, we conclude that for isotopic 
binary mixtures the time scales at which the series in Eq. (1171) can be supposed to be useful 
are the same as those for a single LJ fluid, i.e., of the order of a picosecond. 
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FIG. 3: The critical time scales tq and r* at which the series in Eq. (|17p for the equilibrium single 
component fluid could be supposed to be useful (cf. Sec-IIV] below Eqs. (fT8|) and ([20]) ) as a function 
of the density p for fixed temperature T = 1.0. Note that the physical time scale in picoseconds 
can be calculated by multiplying both r* and tq by the LJ unit time tlj = 2.16 ps. 



C. Nanofluids 

A nanofluid is a binary mixture of LJ fluid particles (A particles) and nanoparticles (B 
particles). For such a mixture, the time scales Tq and and t§ and need not be 
the same. They were here investigated using the same approach as above, but there are 
additional numerical challenges. First of all, for large B particles, the typical relaxation and 
correlation times (say of the particle velocity) grow with increasing R due to the increased 
inertia of the B particle. As a result, it takes longer to equilibrate such a system, and 
one obtains fewer independent data points per time unit. Secondly, since the B particle 
is already quite large, to surround it with a liquid-like fluid of A particles requires a large 
number of A particles. This increase of the number of particles causes a substantial slow 
down of the simulations. To keep down the number of A particles, one takes as few B 
particles as possible. This contributes to a third difficulty, namely, that for the B particles, 
there are fewer particles to average over, leading to poorer statistics. 

Given these difficulties, fewer runs can be performed in a reasonable time for these systems 
and as a result the error bars on the data for the B particles are substantially larger than 
those for the A particles and of the LJ fluids of the previous sections. Nonetheless, we have 
been able to extract estimates for the timescales at which the series in Eq. (11 7p could be 
supposed to be useful also for these systems. 

For the simulations of the nanofluid, two temperature values were taken: a low temper- 
ature T = 1 (corresponding to 122 Kelvin for Argon) and a high temperature T = 3 (366 
Kelvin, chosen to be closer to room temperature). The simulated system contained Nb = 1, 
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it — Z 


It — 4 
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N B = lc A 

r A 
c 6 

c 4 
n B 


283.6 ± 0.4 

-24524 ± 477 

0.036 ± 0.002 
n nfiK -4- n fv*7 


288.4 ± 0.5 

-24865 ± 554 

(29.3 ± 1.7) x 10~ 6 
( i o 4- n k 1 ! v i n — 6 

^ — J_.Z±U.0I X 1U 


296 ±0.7 
-25369 ± 690 
(0.40 ± 0.03) x 10" 6 
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N B = 2 4 

c 6 
c 4 


284.8 ± 0.7 
-24777 ± 788 
0.036 ± 0.002 
_n nvn -L n n^a 

\J . VJ 1 U _L_ VJ . VJtJO 


293.5 ± 0.9 
-25100 ± 751 
(29 ± 2) x 10~ 6 
(— 1 + fil x 10 -6 


308 ± 1 
-26052 ± 1209 
(0.52 ±0.05) x 10- 6 
f— 2 (1+1 Kl vl D^ 9 


iV B = 3 4 

c 6 

c 4 
r B 
c 6 


289.1 ±0.9 
-25500 ± 1039 
0.041 ±0.002 
-0.13 ±0.10 


300 ± 1 
-25717 ± 1170 
(41 ± 2) x 10~ 6 
(-2.2 ±1.2) x 10~ 6 


314 ±2 
-26368 ± 1257 
(0.78 ± 0.06) x 10" 6 
(-3 ±2) x 10~ 9 



TABLE I: The coefficients c\ and Cg for the LJ particles (^4) and the nanoparticles (B) in the 
nanofiuid of Sec. lVI Cl at T = 1. 

2 or 3 nanoparticles of size R = 2, 4 or 6 (i.e., all nine combinations were studied). The 
linear box size was L = 30 so that the number density of the nanoparticles had the values 
p B = 3.7 x lCr 5 , 7.4 x lCr 5 and 1.1 x 1CT 4 for N B = 1, 2 and 3, respectively. To keep the 
properties of the LJ fluid in which the nanoparticles are suspended constant, the remainder 
of the box was filled with LJ particles with a fixed number density pa = Na/ (L 3 — ^7tR 3 N b ), 
which was, somewhat arbitrarily, chosen to be 0.49, i.e. iV^ was chosen such that for given 
L, R and N B , p^ was as close to 0.49 as possible. This required between Na = 11, 912 and 
Na = 13, 227 fluid LJ particles, depending on R and iV^- Note that even though the number 
densities of the nanoparticles are small, by assigning a volume |vri? 3 to each nanoparticle, 
one sees that the volume fraction ranges from 0.124% to 10%. This is a realistic range, 
as experimental volume fractions are of the order of 1%[1]. We did not investigate much 
higher volume fractions to avoid possible complicating effects such as aggregation of the 
nanoparticles. 

For the systems with 1 nanoparticle, 100 runs were performed for each of the two temper- 
ature values T = 1 and T = 3, where first the system was equilibrated using an isokinetic 
Gaussian thermostat, and then the system was run for 8 ps during which the quantities 
appearing in Eqs. fl36|) and fl37j) were measured. For the systems with iV^ = 2, 50 runs 
were performed and for those with N B = 3 the number of runs was 34 (for each temper- 
ature value). Because of the isokinetic Gaussian thermostat, the average over these runs 
approximates the average over the canonical distribution in Eq. 

The resulting values for c\ and Cg are shown in Tables [U and [III for T = 1 and T = 3, 
respectively. From c\ and Eq. (|18p we find the timescales Tq, which are listed in Tables III 1 1 
and IIVI for T = 1 and T = 3, respectively. In Tables HI and HT1 one notices the large error 
estimates for Cg 3 (whose values are negative as in the pure LJ case), which may seem to 
make it hard to draw conclusions from those data. However, according to Eq. (1371) we only 
need the square root of this number to estimate t^, leading to a reduction of the relative 
error by one half, which explains why the results for rf given in Tables II I II and IIVI are still 
reasonable order of magnitude estimates for all cases except for the combination of physical 
parameters R = 6 and T = 3. 
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R = 2 


R = A 


R = 6 


N B = lc A 

c 6 

r B 
c 4 

r B 


9042 ± 7 
(-6.4 ±0.12) x 10 6 
17 ±1 
-442 ± 291 


9034 ± 7 
(-6.3 ±0.11) x 10 6 

(22 ± 2) x 10~ 3 
(-11 ± 17) x 10~ 3 


9123 ±6 
(-6.4 ±0.11) x 10 6 
(357 ± 54) x 10~ 6 
(-18 ±76) x 10~ 6 


N B = 2 cf 

r A 
c 6 

r B 

r B 
C 6 


9049 ± 8 
(-6.4 ±0.13) x 10 6 
17 ±1 
-434 ± 192 


9163 ± 8 
(-6.5 ±0.14) x 10 6 

(23 ± 2) x 10~ 3 
(-13 ± 17) x 10~ 3 


9296 ± 8 
(-6.6 ±0.12) x 10 6 
(390 ± 43) x 10~ 6 
(-18 ±62) x 10~ 6 


iV B = 3 cf 

„A 
c 6 

c 4 
c 6 


9087 ± 8 
(-6.4 ±0.12) x 10 6 
17 ±1 
-425 ± 198 


9203 ± 6 
(-6.5 ±0.1) x 10 6 

(22 ± 1) x 10~ 3 
(-11 ± 11) x 10" 3 


9455 ± 160 
(-6.8 ±0.7) x 10 6 
(468 ± 110) x 10~ 6 
(-4.3 ± 105) x 10~ 6 



TABLE II: The coefficients c 4 and Cg for the LJ particles (A) and the nanoparticles (B) in the 
nanofiuid of Sec. lVI Cl at T = 3. 

We see from Tables ILTT1 and HVl that for the LJ fluid particles (A) surrounding the nanopar- 
ticles, both time scales t a and t a (below which which the expansion of the Van Hove self- 
correlation function around a Gaussian as in Eq. f|T7|) may be useful) are on the order of one 
or two picoseconds. While they decrease moderately with increasing temperatures, these 
time scales are relatively insensitive both to the radius and to the density of the nanoparti- 
cles, and are in fact close to their values in the absence of nanoparticles (cf. Fig.[2]), which 
were also on the order of one to two picoseconds. 

In contrast to this, Tables [TTT1 and HVl shows that the time scales below which the expansion 
of the Van Hove self-correlation function of the nanoparticles (B) around a Gaussian could 
be supposed to be practicable, is considerably larger than for the fluid particles, and, in 
fact, increases with the radius of the nanoparticles up to as much as a factor five for T = 3 
and a factor ten for T = 1 for the largest nanoparticle size studied. The timescales decrease 
upon increasing the density of the nanoparticles, but by a lesser amount, so that the overall 
timescale below which Eq. (1171) is useful is still on the order of five picoseconds for T = 3 
and on the order of ten picoseconds for T = 1 . 

VII. CONCLUSIONS 

We have investigated the short time behaviour of the Van Hove self-correlation function. 
According to Eq. fll7p . for short times, the Hove self- correlation function can be expressed as 
a Gaussian plus corrections, which are proportional to increasing powers of t. For short times, 
this can be re-expressed by the series in Eq. (fTTl) . which is useful provided the contributions 
of the correction terms are small. From the form of these correction terms in Eq. (1171) . one 
sees that they are small at time scales smaller than some critical time scale (tq)- In this 
paper, this time scale was investigated for a number of LJ and LJ-based systems. We found 
that a decrease of the magnitude of the terms in the series Eq. (11 7p occurs below and up to 
the picosecond time scales for LJ fluid particles and up to the ten picosecond time scale for 
nanoparticles. 
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0.763 ± 0.001 
0.833 ± 0.008 
2.40 ±0.03 
i.y ± u.o 
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0.830 ± 0.013 
2.40 ± 0.03 
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0.756 ± 0.001 
0.838 ± 0.013 
5.29 ±0.09 
5 9 + 16 


0.747 ±0.001 
0.842 ±0.019 
7.91 ±0.19 
8 + 3 


N b = ?>t a 
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T G 
T B 


0.759 ± 0.001 
0.825 ±0.017 
2.32 ±0.03 
1.5 ±0.6 


0.752 ±0.001 
0.837 ±0.019 
4.85 ± 0.06 
4.1 ± 1.1 


0.744 ± 0.001 
0.84 ± 0.02 
7.15 ±0.14 
8±3 



TABLE III: The time scales Tq and t£ (in LJ units) as follow from c\ and Cg according to Eqs. (jTHJ) 
and (|20p for the LJ particles (^4) and the nanoparticles (B) in the nanofluid of Sec. lVI Cl at T = 1, 
respectively. Note that the physical time scale in picoseconds can be calculated by multiplying r* 
and tq by the LJ unit time t^j = 2.16 ps. 



Two time scales were in fact calculated: one, denoted by tq, estimates when the first 
correction term to the Gaussian distribution will be small, and the other, denoted by r*, 
estimates the time at which the second correction term is as big as the first one. The 
larger these time scales, the better, since this means that the expansion in Eq. (1171) . i.e., the 
Gaussian plus two correction terms, or perhaps even just the simple Gaussian prefactor, can 
be used for all time scales below (and possibly up to) tq and r*. Note that if these time 
scales are of similar order of magnitude, as they turned out to be, then they could also be 
viewed as a possible estimate of the radius of convergence of the series in Eq. ffTTl) . 

We first investigated the coefficients for the equilibrium pure LJ fluid as a function of 
temperature and concluded that both time scales tq and r* are reduced as a function of in- 
creasing temperature from about 2 picoseconds to 1 picosecond. As a function of density, our 
two estimates of the time scales behave differently. While tq decreases by moderate amounts 
with increasing density, r* stays roughly the same. In all cases though, the timescales are of 
the order of a picosecond or more. One can qualitatively understand the decreasing trend of 
the 'Gaussian' time scale tq for increasing densities, by realizing that the forces between the 
particles perturb the short time ballistic motion away from its Gaussian character. Since 
the forces are stronger at higher densities, the deviations from Gaussian behaviour will then 
occur earlier. 

In mixtures, there is a Van Hove self-correlation function for each component, and cor- 
respondingly, the time scales depend on the component whose Van Hove self-correlation 
function is studied, which is represented by a superscript A = A or B on t g and r*. We 
deduced for a binary isotopic mixture that the time scales Tq and on which Eq. (I17p could 
be supposed to be useful, simply scale as the square root of the mass m\ of the component 
A. As said before, since in Nature, isotopes do not have very large mass ratios, for isotopic 
binary mixtures the time scales at which the series in Eq. f|T7|) is useful are of the same order 
of magnitude as for a one-component fluid, i.e., of the order of a picosecond. 



Finally, we studied these time scales in a recently proposed model of a nanofluid 21 
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0.5560 ± 0.0001 
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0.5561 ± 0.0001 
0.508 ± 0.004 
1.75 ±0.04 
z.4 zb l.y 


0.5547 ±0.0001 
0.508 ± 0.004 
2.7 ±0.1 

A 1 Q 

4 ± o 


N b = 2t£ 

' * 

T G 

T B 

T * 


0.5559 ± 0.0001 
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0.5541 ± 0.0001 
0.504 ±0.005 
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9 9 + 15 


0.5547 ±0.0001 
0.508 ± 0.004 
2.62 ±0.07 
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N B = 3r A 
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0.5553 ± 0.0001 
0.507 ±0.005 
0.89 ±0.01 
0.90 ±0.21 


0.5535 ± 0.0001 
0.504 ± 0.004 
1.74 ±0.02 
2.4 ± 1.2 


0.5498 ± 0.0001 
0.50 ±0.25 
2.5 ±0.15 
9 ±114 



TABLE IV: The time scales Tq and (in LJ units) as follow from c\ and Cg according to Eqs. (jTHJ) 
and (|20p for the LJ particles (A) and the nanoparticles (B), respectively, in the nanofluid of 
Sec. lVI CI at a temperature of T = 3. Note that the physical time scale in picoseconds can be 
calculated by multiplying r* and tq by the LJ unit time tlj = 2.16 ps. 



and found that the time scales are there of the order of five to ten picoseconds for the 
nanoparticles (decreasing with temperature and increasing with radius), while for the fluid 
particles in that model the time scale is still on the order of a picosecond. The difference in 
time scales could be due to the larger mass of the nanoparticles, causing the forces to have 
less influence on their velocities, which therefore remain close to their original (Gaussian) 
distribution for a longer time than in a LJ fluid. It is then no surprise that the distribution 
of displacements for nanoparticles can be described by a Gaussian at longer time scales than 
for the lighter fluid particles. 

One may wonder whether the time scales found in this paper are not so short that 
the classical description on which they were based breaks down. A simple estimate of the 
time scale at which appreciable quantum effects can be expected is given by h/ksT, where 
h is Planck's constant divided by 2ir. At room temperature, this is equal to about 25 
femtoseconds. Note that all of the time scales found in this paper were at picosecond or at 
tens of picosecond scales, i.e., well above this quantum time scale. 

Although our results for the time scales Tq and are only estimates, they are encour- 
aging for the possible application of a Green's function approach to small scale nanometre 
length and picosecond time scales, since the Van Hove self-correlation functions are equilib- 
rium versions of Green's functions!!!, H Furthermore, it is expected that the 
time scales for nonequilibrium systems are similar to those of equilibrium systems, which 
were on the order of picoseconds for fluid particles and on the order of ten picoseconds for 
nanoparticles. This suggests that expansions of the form in Eq. (j!7p can be useful for the 
Green's function approach for transport problems taking place at and below picosecond time 
scales and at nanometre length scales in equilibrium and near-equilibrium systems. 
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APPENDIX A: MOMENTS AND CUMULANTS 

In this appendix we will briefly recall the definitions of the moments and cumulants, and 



how they are related. For more details, see Ref.l24. 

We first remark that multivariate moments and cumulants are simply moments and cu- 
mulants of more than one variable. In general, (multivariate) moments can be defined 
as follows. For a single random variable x with a distribution fi(x), the nth moment is 
«„ = (x n ) = f dx x n fi(x), while for a pair of random variables x\ and x 2 with a joint 
distribution /^(^i, 2^2), the bivariate moments are (x^x^ 2 ) = f dx± f dx 2 x^x^ 2 f 2 (xi, x 2 ), 
and so on for multivariate moments (x™ 1 • • • Xq q ) = f dx± ■ ■ ■ f dx q x™ 1 ■ ■ ■ Xq q fq(x±, . . . , x q ). 
One defines the order of a multivariate moment as the sum Y2l=i n r- F° r near-Gaussian 
(multivariate) distributions, the cumulants are a more convenient way to characterize the 
distribution than the moments, because the cumulants of order higher than two are zero for 
a pure Gaussian. For a single variable the general expression for the nth cumulant K n in 
terms of moments fik<n is 



1 Pe 



—n\ 



oo oo r /m~\'i 

e (E^-O'n^- (AD 

{pe>0} t=l l=\ ^ 

Efci t-Pi=n 

In analogy with the notation u n = (x n ) for moments of a random variable x, one often uses 
the notation n n = ((x n )) for its cumulant sjH]. Here, the superscript n inside the double 
brackets is not a power, as the example ((x 2 )) = (x 2 ) — (x) 2 shows. To avoid confusion, we 
denote instead the cumulants as ((x^)). Therefore, instead of Eq. (lAlj) we can write 



E (f ft -i)in^. <A2) 



{pi>o} e=i 



One can interpret the superscript n between square brackets in this expression as the number 
of 'repetitions' of x. Then, as an alternative to Eq. flA2j) . one can define the cumulants 
recursively as the average of the product of these repetitions minus the product of lower 
order cumulants of all possible groupings of the n repetitions. For instance, for the third 
order cumulant of the displacement one can write ((x' 3 ')) = (x 3 ) — 3((x))((x^)) — ((x)) 3 , 
where the factor three arises from the three ways in which one can group three repetitions 
into a pair and a single repetition. This expression contains the second order cumulant 
((x' 2 l)), which can be written as ((x' 2 ')) = (x 2 ) — ((x)) 2 , while finally ((x)) = (x), leading to 
«4 3 )(t))) = (x 3 }- 3(x)(x 2 } + 2(x} 3 . This is a special case of the general formula flA2j) . 
Similarly to this univariate case, multivariate cumulants can be represented in terms of 
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the averages, in the following way (l8| : 



[x\ x™ 1 )) = -n x \...n, 



{p {e} >o} W W 



(x[\..x q ") \ Pie} 
p {e} \ V " h\...l q \ 



(A3) 



In this notation for the cumulants, quantities separated by semicolons are treated as separate 
random variables and, as above, if a quantity has a superscript within square brackets, it 
denotes that particular number of repetitions of the quantity. Some examples of multi-variate 
cumulants in terms of multi-variate moments are 

«*!» = (sx) (A4) 
({xi,x 2 )) = (xix 2 ) - (xi) (x 2 ) (A5) 
((xi; x 2 \ x 3 )) = (xiX 2 x 3 ) - {x x x 2 ) (x 3 ) - (xix 3 ) (x 2 ) - (xi) (x 2 x 3 ) + 2 (x ± ) (x 2 ) (x 3 ) (A6) 

In the main text, the moments /i and cumulants K occurs as moments and cumulants of 
the displacements of a single particles of a specific component A in a time t, and therefore 
appear with a superscript A (and an implicit time argument t). Furthermore, multi-variate 

cumulants appear where the x 7 are replaced by Y\ 1} or by derivatives of the potential, i.e. 
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